/*
 #
 #  Files       : Matrix.h
 #                ( C++ header file )
 #
 #  Description : The SmallMatrix Library
 #                ( http://code.google.com/p/smallmatrix )
 #
 #  Copyright   : Olivier Juan
 #                ( http://www.mas.ecp.fr/vision/Personnel/juan/ )
 #
 #  License     : CeCILL-C
 #                ( http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html )
 #
 #  This software is governed by the CeCILL-C license under French law and
 #  abiding by the rules of distribution of free software.  You can  use,
 #  modify and or redistribute the software under the terms of the CeCILL-C
 #  license as circulated by CEA, CNRS and INRIA at the following URL
 #  "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and  rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #  same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL-C license and that you accept its terms.
 #
*/

#ifndef _MATRIX_H
#define _MATRIX_H

#include <SmallMatrix/SmallMatrix.h>

NAMESPACE_BEGIN(SmallMatrix)

template <int M,int N> class Matrix: public Matrix_Base<Matrix<M,N>,M*N,M,N>
{
public:
	Matrix():Matrix_Base<Matrix<M,N>,M*N,M,N>(){};
	Matrix(const reel &v):Matrix_Base<Matrix<M,N>,M*N,M,N>(v){}
	Matrix(const reel tab[M*N]):Matrix_Base<Matrix<M,N>,M*N,M,N>(tab){}
	template <typename TT,int SSIZE> Matrix(const Matrix_Base<TT,SSIZE,M,N>& copy):Matrix_Base<Matrix<M,N>,M*N,M,N>(copy){}
	inline reel& operator()(const int& idxm,const int& idxn);
	inline const reel& operator()(const int& idxm,const int& idxn)const;
	template <typename T1,int SIZE1,int MN,typename T2,int SIZE2> inline Matrix<M,N>& prod(const Matrix_Base<T1,SIZE1,M,MN>& M1,const Matrix_Base<T2,SIZE2,MN,N>& M2);
	template <typename T1,int SIZE1,typename T2,int SIZE2> inline Matrix<M,N>& add(const Matrix_Base<T1,SIZE1,M,N>& M1,const Matrix_Base<T2,SIZE2,M,N>& M2);
	template <typename T1,int SIZE1,typename T2,int SIZE2> inline Matrix<M,N>& subst(const Matrix_Base<T1,SIZE1,M,N>& M1,const Matrix_Base<T2,SIZE2,M,N>& M2);
	template <typename T1,int SIZE1> inline Matrix<M,N>& add(const Matrix_Base<T1,SIZE1,M,N>& M1);
	template <typename T1,int SIZE1> inline Matrix<M,N>& subst(const Matrix_Base<T1,SIZE1,M,N>& M1);
	template <typename TT,int SSIZE> inline Matrix<M,N>& mult(const Matrix_Base<TT,SSIZE,M,N>& Two);
	template <typename TT,int SSIZE> inline Matrix<M,N>& div(const Matrix_Base<TT,SSIZE,M,N>& Two);
	Matrix<N,M> trans()const;
};
//! Read/Write Accessor
template <int M,int N> inline reel& Matrix<M,N>::operator()(const int& idxm,const int& idxn)
{
	assert((idxm>=0)&&(idxm<M)&&(idxn>=0)&&(idxn<N));
	Matrix<M,N>& ref=*this;
	return ref[idxm*N+idxn];
}
//! Read Only Accessor
template <int M,int N> inline const reel& Matrix<M,N>::operator()(const int& idxm,const int& idxn) const
{
	assert((idxm>=0)&&(idxm<M)&&(idxn>=0)&&(idxn<N));
	const Matrix<M,N>& ref=*this;
	return ref[idxm*N+idxn];
}
//! Matrix-Matrix Product in which the result is a Matrix
template <int M,int N> template <typename T1,int SIZE1,int MN,typename T2,int SIZE2> inline Matrix<M,N>& Matrix<M,N>::prod(const Matrix_Base<T1,SIZE1,M,MN>& M1,const Matrix_Base<T2,SIZE2,MN,N>& M2) {
	Matrix<M,N>& ref=*this;
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++) {
			ref(i,j)=0;
			for (int k=0;k<MN;k++)
				ref(i,j)+=M1(i,k)*M2(k,j);
		}
	return ref;
}
//! Matrix-Matrix Sum in which the result is a Matrix
template <int M,int N> template <typename T1,int SIZE1,typename T2,int SIZE2> inline Matrix<M,N>& Matrix<M,N>::add(const Matrix_Base<T1,SIZE1,M,N>& M1,const Matrix_Base<T2,SIZE2,M,N>& M2) {
	Matrix<M,N>& ref=*this;
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++)
			ref(i,j)=M1(i,j)+M2(i,j);
	return ref;
}
//! Matrix-Matrix Substract in which the result is a Matrix
template <int M,int N> template <typename T1,int SIZE1,typename T2,int SIZE2> inline Matrix<M,N>& Matrix<M,N>::subst(const Matrix_Base<T1,SIZE1,M,N>& M1,const Matrix_Base<T2,SIZE2,M,N>& M2) {
	Matrix<M,N>& ref=*this;
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++)
			ref(i,j)=M1(i,j)-M2(i,j);
	return ref;
}
//! Matrix-Matrix Sum in which the result is a Matrix : add M1 to the Matrix
template <int M,int N> template <typename T1,int SIZE1> inline Matrix<M,N>& Matrix<M,N>::add(const Matrix_Base<T1,SIZE1,M,N>& M1) {
	Matrix<M,N>& ref=*this;
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++)
			ref(i,j)+=M1(i,j);
	return ref;
}
//! Matrix-Matrix Substract in which the result is a Matrix : substract M1 to the Matrix
template <int M,int N> template <typename T1,int SIZE1> inline Matrix<M,N>& Matrix<M,N>::subst(const Matrix_Base<T1,SIZE1,M,N>& M1) {
	Matrix<M,N>& ref=*this;
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++)
			ref(i,j)-=M1(i,j);
	return ref;
}
template <int M,int N> template <typename TT,int SSIZE> inline Matrix<M,N>& Matrix<M,N>::mult(const Matrix_Base<TT,SSIZE,M,N>& Two) {
	Matrix<M,N>& ref=*this;
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++)
			ref(i,j)*=Two(i,j);
	return ref;
}
template <int M,int N> template <typename TT,int SSIZE> inline Matrix<M,N>& Matrix<M,N>::div(const Matrix_Base<TT,SSIZE,M,N>& Two) {
	Matrix<M,N>& ref=*this;
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++) {
			assert(Two(i,j)!=0);
			ref(i,j)/=Two(i,j);
		}
	return ref;
}
//! Matrix Transpose
template <int M,int N> Matrix<N,M> Matrix<M,N>::trans()const {
	const Matrix<M,N>& ref=*this;
	Matrix<N,M> res;
	for (int i=0;i<N;i++)
		for (int j=0;j<M;j++)
			res(i,j)=ref(j,i);
	return res;
}


NAMESPACE_END

#endif
